% % Obtain statistics of simulations of the two-productivity model in Section 5
% Octoer 2022
% Makoto Nirei

%load simu.mat;
dlogPm=zeros(maxmonth,pilength); L=zeros(maxmonth,pilength);    % monthly aggregate price increase
dlogPm(1,:)=sum(dlogP(1:calvom(1),:));                  % monthly price increase for the 1st month
idl=(nump(1:calvom(1),:)>0);                            % indicator for price rise
idlm=ones(calvom(1),pilength)-idl;                      % indicator for price cuts
cumcost_t = cumsum(nump.*(nump>1)*delta - nump.*(nump<-1)*delta_u); % cumulative menu costs
DL_t = cumcost_t -interp1([0; calvot],[zeros(1,pilength); cumcost_t],calvot-1/12,'previous'); % cumulated in the last month. 1st month will be discarded
mNpY=zeros(maxmonth-1,pilength);
mcons=mNpY; DL=mNpY;
pc1=(1-alpha_U*(1-gamma_U))/(1-alpha_U)/(1-gamma_U);
m_ss=zeros(1,pilength); c_ss=m_ss; n_ss=m_ss; w_ss=m_ss;
for i_pi=1:pilength
    ss=ss_pi{i_pi};
    m_ss(1,i_pi)=ss.m;
    c_ss(1,i_pi)=ss.C;
    n_ss(1,i_pi)=ss.N;
    w_ss(1,i_pi)=ss.w;
end

options = optimset('Display','off');
for i=2:maxmonth
    dlogPm(i,:)=sum(dlogP(calvom(i-1)+1:calvom(i), :)); % monthly price increase
    idl=(nump(calvom(i-1)+1:calvom(i),:)>0);            % indicator for price rise
    idlm=ones(calvom(i)-calvom(i-1),pilength)-idl;      % indicator for price cut
    DL(i-1,:)=DL_t(calvom(i-1)+1,:);
    mNpY(i-1,:) = (sum(NpY(calvom(i-1)+1:calvom(i),:).*diff(calvot(calvom(i-1)+1:calvom(i)+1)))...
        +NpY(calvom(i-1),:)*(calvot(calvom(i-1)+1)-month_time*(i-1))...
        -NpY(calvom(i),:)*(calvot(calvom(i)+1)-month_time*i))/month_time; % monthly averaged NpY
    cc0=(1-DL(i-1,:)/n/month_time)./mNpY(i-1,:);
    cc1=cc0./(m_ss*rho*alpha_U).^(1/(1-alpha_U)/(1-gamma_U));
    compu_c = @(x) x -cc0 + cc1.*x.^pc1;
    mcons(i-1,:)=fsolve(compu_c, c_ss, options);
end

dlogPy=zeros(maxmonth-11,pilength);
for i=1:maxmonth-11
    dlogPy(i,:)=sum(dlogPm(i:i+11,:));                  % monthly, year-on-year price increase
end
stdP = std(dlogPy);                                     % volatility of annual inflation
meanP = mean(dlogPy);                                   % average of annual inflation

% stats for various pi
disp('pi, mean(dP/dt), std(dP/dt), downward repricings');
disp([piv; meanP; stdP; negnump]);

disp('std of log(c) for various pi');
disp(std(log(mcons)));

disp('std of monthly log(N/Y)');
disp(std(log(mNpY)));

disp('mean, std, min, max of monthly d(log(P))');
disp([mean(dlogPm(:,:)); std(dlogPm(:,:)); min(dlogPm(:,:)); max(dlogPm(:,:))]);

i_pi=2;
figure(11),
set(gca,'FontSize',20);
stairs(calvot(100:200)-calvot(100),cumsum(dlogP(100:200,i_pi))+1);
xlabel('Time (unit time is a year)');
ylabel('log P');
xlim tight;
title('Typical path of log P(t)');

figure(12),
set(gca,'FontSize',20);
plot(0:100,dlogPm(100:200,i_pi));
xlabel('Months');
ylabel('logP(t+1)-logP(t)');
title('Typical path of monthly inflation');

disp('Autocorrelations for dlogP and dlogPm');
corrcoef(dlogP(100:end-1,2),dlogP(101:end,2))
corrcoef(dlogPm(100:end-1,2),dlogPm(101:end,2))

disp('statistics used in Appendix D');

mN(1:maxmonth-1,:)=mNpY(1:maxmonth-1,:) .* mcons(1:maxmonth-1,:)./(1-DL(1:maxmonth-1,:)/n/month_time);
mwage(1:maxmonth-1,:)=(1-alpha_U)/alpha_U *mcons(1:maxmonth-1,:)./(1-mN(1:maxmonth-1,:));

disp('stdev and autocorr of log(w)');
std(log(mwage(:,:)))
corrcoef(mwage(1:end-1,i_pi),mwage(2:end,i_pi))
[h,pValue,stat,cValue]=adftest(mwage(:,2));

disp('V(L)/E(L)');
Lpos=nump(nump(:,2)>1, 2)-1; % subtracting the Calvo
var(Lpos)/mean(Lpos)